{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sisl\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this example, you will familiarize your-self to the concept of *buffer atoms*. A buffer atom is an atom that is completely neglected in the TranSiesta self-consistent calculation, but is used as an initialization for the bulk electrode regions. \n",
    "\n",
    "Here a pristine graphene flake will be constructed, and subsequently a Carbon chain will act as an STM-like tip to simulate STM experiments.\n",
    "\n",
    "As the Carbon chain is terminated to vacuum, the dangling bonds will create spurious effects very different from a pristine, bulk chain. To understand why it is necessary to add buffer atoms, it is useful to understand the TranSiesta method. **Any** TranSiesta calculation starts with calculating an initial guess for the Hamiltonian as input for the Green function method:\n",
    "\\begin{equation}\n",
    "   \\mathbf G^{-1}(E) = \\mathbf S (E+i\\eta) - \\mathbf H - \\sum_i\\boldsymbol\\Sigma_i\n",
    "\\end{equation}\n",
    "If the initial $\\mathbf H'$ represents a Hamiltonian close to the open-boundary problem $\\mathbf H$; it will converge with a higher probability, and in much less time. Improving the initial guess Hamiltonian is time well-worth spent as TranSiesta is, typically, more difficult to converge. The initial guess Hamiltonian is a Siesta calculation with full periodicity.\n",
    "\n",
    "As an example, consider the Hamiltonian for the chain:\n",
    "\n",
    "     <vacuum>     C -- C -- C -- C -- C -- C ...\n",
    "     \n",
    "It is clear that the atom closest to the vacuum region resides in a *very* different chemical and potential landscape than an atom in the middle of the chain. If TranSiesta uses the initial Hamiltonian for the chain electrode as the atom closest to the vacuum region, it will be very far from the potential landscape of a bulk electrode. So to mitigate this one can specify:\n",
    "\n",
    "    %block TBT.Atoms.Buffer\n",
    "      atom [ 1 -- 2 ]\n",
    "    %endblock\n",
    "\n",
    "to remove the first 2 atoms from the TranSiesta calculation (note that negative indices counts from the *end*). Then the electrode will begin from the 3rd atom, which is farther from the dangling bond. This will be a *much* better initial guess for the Hamiltonian. Other strategies to improve the potential landscape may be to terminate the dangling bonds with Hydrogens or other atomic species.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "graphene = sisl.geom.graphene(1.44)\n",
    "elec = graphene.tile(2, axis=0)\n",
    "elec.write('ELEC_GRAPHENE.fdf')\n",
    "elec.write('ELEC_GRAPHENE.xyz')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "C1d = sisl.Geometry([[0,0,0]], graphene.atoms[0], [10, 10, 1.4])\n",
    "elec_chain = C1d.tile(4, axis=2)\n",
    "elec_chain.write('ELEC_CHAIN.fdf')\n",
    "elec_chain.write('ELEC_CHAIN.xyz')\n",
    "chain = elec_chain.tile(3, axis=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "device = elec.tile(5, axis=1).tile(4, axis=0)\n",
    "\n",
    "# Attach the chain on-top of an atom\n",
    "# First find an atom in the middle of the device\n",
    "idx = device.close(device.center(what='xyz'), R=1.45)[1]\n",
    "# Attach the chain at a distance of 2.25 along the third lattice vector\n",
    "device = device.attach(idx, chain, 0, dist=2.25, axis=2)\n",
    "# Add vacuum along chain, we really no not care how much vacuum, but it\n",
    "# is costly on memory, not so much on performance.\n",
    "device = device.add_vacuum(15, axis=2)\n",
    "device.write('DEVICE.fdf')\n",
    "device.write('DEVICE.xyz')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "- Add missing electrode information in `RUN.fdf`.\n",
    "- Perform all required TranSiesta calculations, first the electrodes, then the device region.\n",
    "- Create a new directory for a different range of buffer atoms, from 0 to 4, start by using 4 buffer atoms.  \n",
    "  How does convergence behave for different number of buffer atoms?\n",
    "  \n",
    "  **REMARK** there are 2 places in the fdf file you should change when changing the number of atoms (the electrode atom specification *and* the buffer atoms).\n",
    "- **TIME**: one can combine electrode options `bulk` and `DM-init` to improve the initial $\\mathbf H$ for TranSiesta. Take a system with 1 buffer atom and play with the effect of these options.\n",
    "- Calculate transport properties for all (converged) TranSiesta calculations\n",
    "- Plot the transmission and DOS for all TBtrans calculations, do they differ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Adapt to read in the siesta.TBT.nc from different directories and plot them.\n",
    "tbt = sisl.get_sile('siesta.TBT.nc')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
